'''
Created on Oct 15, 2009

@author: mkiyer
'''

from optparse import OptionParser
from veggie.app.chipseq.macs import MACS
from veggie.app.chipseq.hpeak import HPeak
import logging
import pysam

def run_macs(t_file, c_file, format="SAM", name="macs"):
    # get tag size
    samfile = pysam.Samfile(t_file)
    for read in samfile.fetch():
        read_length = read.rlen
        break
    samfile.close()
    # default options
    macs_options = {'--mfold': 32,
                    '--tsize': read_length, 
                    '--pvalue': '1e-5',
                    '--wig': None,
                    '--verbose': 2}    
    macs_app = MACS(options=macs_options)    
    macs_app.run(t_file, c_file, format, name=name)

def run_hpeak(t_file, c_file, format="BED", name="hpeak"):
    hpeak_options = {'-fmin': 175,
                     '-fmax': 225,
                     '-sig': '1e-5',
                     '-wig': None,
                     '-ann': None}
    hpeak_app = HPeak(options=hpeak_options)
    hpeak_app.run(t_file, c_file, format=format, name=name)

if __name__ == '__main__':
    optionparser = OptionParser("usage: %prog [options]")
    optionparser.add_option("-t", "--treatment", dest="treatment")
    optionparser.add_option("-c", "--control", dest="control", default=None)
    optionparser.add_option("-n", "--name", dest="name")
    optionparser.add_option("-p", "--peakfinder", dest="peakfinder",
                            default="macs")
    (options, args) = optionparser.parse_args()
    
    if options.peakfinder == "macs":
        run_macs(options.treatment,
                 options.control,
                 format="SAM",
                 name=options.name)
    elif options.peakfinder == "hpeak":
        run_hpeak(options.treatment,
                  options.control,
                  format="BOWTIE",
                  name=options.name)
    else:
        print "peakfinder %s not recognized" % options.peakfinder
    